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THE  APPLICATIONS  OF  FINITE-ELEMENT  METHOD 
IN  NUMERICAL  WEATHER  FORECAST 

Mathematics  Department 
Shantung  University 

and 

Shen  ChAng  Sze  and  Chen  Shih  Chieh 
Weather  Bureau,  Shantung  Province 

(received  December  15,  1975) 

The  finite-element  method  has  wide  applications  in  structural  mechanics 

<i,*> 

and  elasticity  mechanics.  In  the  last  two  or  three  years,  the 

possibility  of  applying  finite-element  methods  to  numerical  weather 
forecast  has  been  investigated  from  a theoretical  point  of  view.^' ^ 

This  article  is  a discussion  of  such  applications  based  on  a positive- 
pressure  vorticity  equation.  The  grid  division  used  consists  of  large 
and  small  rectangles.  Such  grids  have  manifested  effectiveness  in  im- 
proving the  forecast  accuracy  of  the  region  under  consideration,  and  they 
are  also  economical  in  terms  of  computer  storage.  Although  rectangular 
grids  are  employed  in  some  foreign  nations,  the  finite-difference  method 
is  always  used  in  the  computation.  In  the  finite-difference  method,  cal- 
culations are  performed  twice  (or  three  times),  namely,  calculations  are 
done  first  on  the  coarse  grid  and  then  on  the  fine  division  grids  locally. 

In  so  doing,  artificial  boundary  conditions  must  be  imposed  on  the 
boundaries  of  the  local  fine  grids.  Such  difficulties  can  be  overcome  by 
the  finite-element  method.  Judging  from  the  calculations  for  the  situations 
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studied,  the  application  of  the  finite-element  method  in  numerical  weather 

forecast  is  a meaningful  undertaking  from  both  practical  and  theoretical 

standpoints.  It  should  be  pointed  out  that  the  finite-element  method 

C5) 

allows  one  to  generalize  the  Arakawa  scheme  , a scheme  usually  used 
in  overcoming  the  nonlinear  instability,  to  the  general  grid  division 
and  obtain  a neat  and  rigorous  proof  of  the  conservation  laws.  Since  the 
energy  conservation  problem  for  an  arbitrary  grid  is  solved,  it  would  be 
eventually  possible  to  carry  out  numerical  forecast  using  monitor  station 
values  directly. 

(4> 

In  the  course  of  our  numerical  calculation,  Jespersen  s work  came 
to  our  attention.  He  also  discovered  that  the  Arakawa  scheme  is  an  Inevitable 
result  of  the  finite-element  method.  However,  Jespersen  considered  only 
the  situation  of  equal  distance  grids.  Hence,  Jespersen's  work  is  a 
special  case  of  our  development. 

I.  Finite-Element  Equation  * 

For  the  sake  of  simplicity,  consider  the  initial  boundary  value 
problem  of  a positive  pressure  vorticity  equation. 

- a C ^f)  + # |f  **+#).  ^ — CO 

HU10  = (a> 

The  notations  employed  here  are  the  ones  commonly  used  in  numerical  weather 
forecast.  G is  the  rectangular  forecast  region  and~^G  is  the  boundary 
of  G. 

This  article  was  received  on  the  5th  of  Decembex  1975. 

^Computational  mathematics  major  grade  72  students  or  the  Mathematics 

Department,  Shantung  University,  participated  in  numerical  calculations. 
Comrade  Chang  Ching  Hua  of  the  Shantung  Province  Weather  Bureau  participated 
in  some  data  reduction  and  graphing. 
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In  order  to  arrive  at  the  finite-element  equations  of  (1)  and  (2), 

(1),  (2) 

they  are  transformed  to  their  equivalent  variation  problems^  That  is, 


among  the  proper  smooth  functions  satisfying  the  boundary  condition 
U J ^ “ 0,  the  general  function  F (II)  given  below  is  minimized  to  the 


function!/  (X,  J)  , where  U 

f(u)  =If{  (fJ^T  +#  f f*  +f) J 


(4) 


Fig.  1 

According  to  the  finite-element  method,  region  G is  divided  into 
a number  of  rectangular  elements.  We  assume  U to  be  continuous  on  G 
and  to  be  a bilinear  function  of  each  rectangle.  The  apices  of  the 
rectangular  element  O e are  i,  j,  m,  and  p,  and  the  lengths  of  the  two 
sides  arc  a and  b (see  Fig.  1).  The  coordinates  of  points  i,  j,  m and 
p are  (XJ,  ft),  (Xj,  Jfj),  (XM,  , and  (X^,  ^ ) respectively.  The 
trans  format ion 

(5) 

S = + bTl " 

changes  the  unit  rectangle  to  a unit  square  in  the  ( S -7  ) plane  with 

(0,  0)  and  (1,  1)  being  its  diagonal  apices.  The  bilinear  interpolation 
function  on  this  element  Q e can  then  be  shown  to  be 


*(.*,  y)  “ (1  — { — n + {i 

+ (i  — {»»)“,  + tn"m  + (n  — {»»)««„ 


(6) 


IWMHMi 


where  U,’,  Uj,  and  Up  are  the  values  of  U (*.  </ ) at  node  points  i,  j, 
m,  and  p,  respectively.  On  Q e,  let 

f « -—4$  + f>  (7) 

functions  Q and  £ can  also  be  expressed  in  terms  of  bilinear  interpolation 
function  similar  to  (6): 

here  subscipted  quantities  are  values  of  the  respective  functions  at  points 
indicated  by  the  subscripts.  One  can  then  write  down: 

T(*,i)sdxSif  5* 

= £ lo-pdi  oi-vcii  v $mo 

(9) 

* o-t- tX*j  rr 

The  set  of  algebraic  equations  satisfied  by  the  values  of  |(  at  the  nodes 
can  be  obtained  by  substituting  (6)  and  (9)  into  (4),  evaluating  the  integral, 
and  then,  according  to  variation  principles,  finding  the  minima  of  all  Ufc. 

This  set  of  equations  is  commonly  known  as  the  finite-element  equation. 


Fig.  2 
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ft 


If  point  0 is  an  inner  node,  the  only  relevant  values  of  u in 
finding  the  equations  with  respect  to  point  0 using  the  method 
described  above  are  those  at  the  nine  corners  of  the  four 
adjacent  rectangles  sharing  point  0 as  their  common  apex.  If 
the  nodes  are  indexed  as  shown  in  Fig.  2,  the  finite-element 
equation  obtained  with  respect  to  point  0 is: 

2*1,1..- *1*1,  (10) 

* •• 

where  the  coefficients  are  listed  below: 


*101  ■ 

*111  ■ 

*121  ■ 

*131  ■ 

*H1  ■ 

*1^1  ' 
*161  < 
*171- 
*l«1  • 
T, 
T, 


+ Ji.  + A.  + A. 

3#,  3*.  is,  34, 

+ A + A.  + A.  + #1 


3«j 


34, 


3«,  n, 


+ -1(T,  + T,+  T.+  T.). 
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A_ 

3*i 

1 


a | 
64, 


A 

3a, 


64, 


+ -CT.  + r,,. 


A 

6a, 


a, 

34, 


A. 

6a, 


A. 

34, 


+ f,(r'  + r')’ 


-A  + -> 


- A.+  • 


3a,  64,  3a, 

+ Jj(r,  + T,). 


OD 

L 

64) 


Si.  + A a- 


_ 

6a,  34j  6a, 

+ ^(T,+  TO, 

_ A__a_  + Ii 

6a,  64,  36’ 


34, 


ii. 

64, 


- *■ 

6a, 

- A--A.+ 

- A.-A.+ 


li 
36’ 

li 
36’ 

li 

64,  36’ 


6a, 

«i4„  T,  — a, 4,. 

m*  ml 

^•i*n  T,  — ^a,4,. 


1 


6 


The  right  hand  side  term  of  Eq.  (10)  is 

b * G0^  + G 0X+)  , (12) 

"he"  $.**-»-  w,  - 0 

| <?;***,(&- s,)-  *,( ?i- -jj+vs-o 

] (13) 

l C -- ?»(#r-4t)-r.(#,-«,)-3:,a-%)+rJ(^.4r) 

If  Eq.  (10)  is  multiplied  by  4/  ^(a^  + a£)  (b^  + b2)J  » it  can  easily 
be  seen  that,  when  grids  are  equal  distant  (a-^  = a2  ■ b^  = b£  = h) , the 
above  equation  is  similar  to  the  Helmholz  equation  with  an  average  nine- 

point  difference  scheme;  furthermore,  the  right  hand  side  term  is  exactly 

(S) 

the  Arakawa  conservation  expression. 


Fig.  3 Large  grid  distance  is  h and  small  grid  distance  is  h/2.  (h  * 540  km) 


We  divide  the  forecast  region  G into  grids  of  intervals  h and  h/2  (h  is 

540  km)  as  shown  in  Fig.  3.  The  finite-element  equation  obtained  with 

the  method  described  above  is  then  solved  to  get  the  values  of  the  function 

at  grid  nodes.  The  conventional  three-step  method  is  then 

used  to  extrapolate  the  forecast  field  ^ at  24  and  48  hours  with  a time 
At 

interval^of  *s  hour.  In  calculating  the  term  on  the  right  hand  side,  the 
values  of  A§  at  the  nodes  are  arrived  at  by  Spline  extrapolation  and 
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the  values  of  ^ at  the  nodes  are  then  obtained.**  The 

results  of  our  numerical  computation  indicated  that  the  grids  are  stable 
and  that  no  spatial  smoothing  was  necessary  during  the  computation.  The 
resultant  48  hour  forecast  diagram  was  relatively  close  to  the  true 
flituatioa.  Contrary  to  the  finite  difference  method  where  additional 
boundary  conditions  are  imposed,  the  finite-element  method  gives  rise  to 
the  equations  at  the  boundary  nodes. 


II.  Generalization  of  the  Arakawa  Conservation  and  Energy  Conservation 


Definition: 


2 !j '7C4.T)Nt‘<lMy 
Z If  Ne' 


(14) 


e'  ae' 

where  De#  is  any  element  containing  point  0 as  an  apex  and  the 
summation  is  over  all  such  elements.  The  factor  Ne#  has  the  following 

convention:  If  Q e<  is  a rectangle  and  point  0 assumes  the  "i"  corner 

of  □ e*  , then  Ne*  = = ( ) ; when  point  0 assumes  the 

"j"  corner  in  element  Q g/  , Ne*  « N£  * ( ^ ^ )»  when  point 

0 is  the  "m"  corner  of  Q g/ , Ne#  » - I1?  , and  finally  when  point 

0 is  the  "p"  corner  of  Q e/  , Ne<  “ N4  * ^ . 


**  The  value  of  p at  the  nodes  can  be  obtained  by  the  following 
equation  if  Spline  extrapolation  value  is  not  used: 

fa)  -4-1  *0(  !>.+  <>»)  1 , 

t** 


A similar  equation  can  be  derived  for  the  boundary  points. 
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Eq.  (14)  is  the  dispersion  formula  of  J(  ^ J ).  The  denominator 
of  Eq.  (14)  as  S^.  If  point  0 is  the  inner  node  shown  in  Fig.  2,  Eq.  (14) 


gives 


j - * + G» * + W* 

3(«,  + «,)(4i  + 4,)’ 


. n ++  +X  _ X+  , - 

where  Gq  , Uq  and  Gq  are  given  in  Eq.  (13).  When  ^ = 4/  L(aj  + 

(b^  + b2)J  , Eq.  (15)  reduces  to  b[o]  / S,  . Thus,  expression  (15)  is 

* W which  is  Eq.  (12) 

consistent  with  the  right  hand  term  of  Eq.  0.0)^ divided  by  and,  obviously, 
it  is  the  generalization  of  Arakawa's  conservation.  When  a^  = a2  = b-^  = b2  = 
h,  Eq.  (15)  becomes  the  usual  Arakawa  conservation  expression: 


+ jr  + jr). 


where 


.♦*  . ..  1 j+»  _ _1_  r+*  /*+«= 

J o “ttt0*  » iUu'  *7° 


Eq.  (15) 


is  the  generalized  Arakawa  conservation  for  nonuniform  rectangular  grids 
and  Eq.  (14)  is  the  generalized  Arakawa  conservation  for  an  arbitrary 


Eq.  (14)  also  satisfied  the  following  three  conservation  relations: 


J-  So  “ 0 (Conservation  of  mean  vorticity) 


Z $J0  S„  = 0 (Conservation  of  mean  kinetic  energy)  (17) 

0 

Z 5-#J0  S„  =0  (Conservation  of  mean  square  vorticity)  (18) 

where  0 assumes  all  nodes  including  nodes  on  the  boundary  and  S.  is  the 

o 

denominator  of  Eq.  (14). 

As  a matter  of  fact,  since  Jocobian  satisfies  the  conservation  of 


mean  vorticity: 


J ($»5  ) dxdy  - Oj 


i 
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we  have 


j(#.r  > 


dxdy  = 0. 


(19) 


Notice  that 


N.  = 1 » 


i = 1 


A 4 

hence,  Jj  J(  4 , "£  ) dxdy  - El!  ■'*  , £ ) N±  dxdy.  (20) 


□e 


i = 1 


From  Eq.  (14),  we  have 


If  J < ^ ’ 5 > N ' 

o O e'  tv 


dxdy. 


and  the  above  expression  can  be  rewritten  by  changing  the  summation  over 


nodes  0 to  a summation  over  all  elements  ft 


2 3e  s,  = Z 2 II j ( <% . r > Ni 


dxdy. 


« t’~'  a 


and  from  (19)  and  (20),  we  have 

^ S0  =0  (Conservation  of  mean  vorticity) . 

o 

The  proof  of  (17)  and  (18)  are  the  same  so  only  the  proof  of  (17) 
will  be  given  below.  From  the  integral  form  expression  of  the  conservation 
of  mean  kinetic  energy,  we  have 


(21) 


= 2.  Jj 

e oc 


= o. 


MW,.  _ 
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Further,  from  (14)  we  have 

LiT.s.'ZiZti  JUrtNe'Mf 

“T^  0 J»  . e’  ot- 

= Z.Z  JT  $.30,'*') N*' 

o «'  oe» 


where  is  the  summation  over  all  rectangles  Q , having  point  0 

e'  6 


as  an  apex.  Again  change  the  summation  over  nodes  to  a summation  over 
each  element  Q e,  it  can  be  seen  that  the  equation  above  becomes: 

Z2>T6S.  =2  I // 

o 6 e Qe 

where  ^ (i  = 1,  2,  3,  4)  represents  the  value  of  ^ at  the  four 


corners  i,  j,  m,  and  p of  the  element  Q e:  - £i,  <2  = £ j. 


$ = £ and  & = £ . 

J-  m ■*  e I,  o 


e 2 


3 U1  'A 

it  is  evident  that  on  the  element  £3  > we  have 


From  the  bilinear  interpolation  equation  (8), 


4 

1*1 


Ni  " "1  + "2  + N3  + "4  • £ 


Substituting  the  above  expression  into  (22)  and  using  Eq.  (21),  we  obtain 

Sc  = Z il  £ J ( 5 . 5 > dxdy  = 0 

° * 

and  thus  completed  the  proof  of  (17).  Eq.  (18)  can  be  proven  in  a similar 
manner. 

In  conclusion,  the  Jc  defined  in  Eq.  (14)  satisfies  relations  (16) 
through  (18).  Generalization  of  this  kind  cannot  be  achieved  by  the 
difference  method;  furthermore,  the  proof  of  the  three  conservation 
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relationships  using  the  finite-element  method  is  simpler  and  more  rigorous 
as  compared  to  the  proof  of  Arakawa  even  for  the  case  of  equal  distance 
grids  where  conservations  can  be  proven  by  the  difference  method. 
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